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Abstract 


In this thesis, various aspects of modelling of non-st ationary signals are considered. 

A variety of non-parametric and parametric methods has been introduced in lit- 
erature to analyse non-stationary signals. The present thesis deals with relating 
various models to a unified framework and addressing key issues of non-stationary 
modelling. This, in turn, consolidates the existing techniques and finds new avenues 
in parametric modelling of non-stationary signals. 

The non-stationary signals that we come across in nature vary over time in the 
following aspects, viz., in frequency content or in energy content or in both. For 
estimation of evolutionary spectrum by parametric rational transfer function model, 
so far, attention was focused only on the time variation of AR/ARMA coefficients. 
The time variation of gain (i.e. the input noise energy) was ignored in this regard. 
In the present work, it is illustrated that both of those features are important in 
their own ways. The time variation of AR co-efficients indicates the change in 
spectral shape with time, and the time-varying gain stores the information about 
signal intensity. Moreover, the time variation of one feature can not compensate for 
the time variation of other feature in the estimated evolutionary spectrum. While 
modelling non-stationary speech signals by time-dependent AR model, the time 
variation of gain becomes prominent. This feature is really a characteristic of the 
production mechanism of speech. We present here a non-iterative technique to 
estimate the gain parameter. Our study reveals that the AR model with time- 
varying gain is a suitable model for speech phoneme. 

Various types of parametric models for non-stationary signals are in use at 
present. These are the damped sinusoid or complex exponential model, amplitude 
modulation (AM), frequency modulation (FM) and amplitude/frequency modula- 
tion (AFM) models. Damped sinusoids are used for modelling transient signals like 
nuclear magnetic resonance (NMR) and electromagnetic pulse (EMP) responses. It 
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is known that the speech signal can be generated by the AFM model. In radar and 
sonar applications, the FM and AFM signals are encountered. 

The main diffi culty that we face in employing the above mentioned modulated 
processes (AM, FM and AFM) to model non-stationary signals is that it leads 
to a set of non-linear equations, which are to be solved. Our study reveals that 
all the above mentioned models are only the special cases of the time-dependent 
ARMA (TARMA) model. Remember here how an undamped sinusoidal signal can 
be modelled by a constant coefficient AR process and the signal parameters can be 
estimated through an estimation of the process parameters. It is shown that in the 
general fr am ework these models keep their own distinctive features, and they can 
be class if ied according to the time dependence of position/locus of the pole of the 
TARMA model. Since various models can be related to the general TARMA model, 
the equivalence provides a unified approach for synthesis and analysis of the signal. 
By using a set of basis functions to represent the time variations of the parameters 
of the TARMA model, one may transform the convoluted non-linear estimation 
problems to the simple linear ones. This new approach of parametric estimation 
increases the accuracy of estimation of the base-band signal, compared to what can 
be achieved by the existing non-parametric methods. The introduction of time- 
dependent damping factor adds a new dimension in modelling transient data in 
terms of accuracy /flexibility of the model, as shown in the thesis. 

The ECG is an example of a non-stationary signal with its characteristic pseudo- 
periodicity. The amplitude modulated sinusoidal model is a natural choice for such 
signal. It is found that when the modulating signal is exponential function of 
time, it can represent the burst type QRS complex very well. This model not only 
can synthesise the ECG signal as a whole but also its constituent waves. As the 
constituent waves of ECG give indications of various activities of the heart, this 
model may be useful in diagnosis. 

The thesis is written to highlight the need of modelling any natural signal accord- 
ing to its own characteristics. The study shows effectiveness and area of application 
for each modelling technique, and how estimation of parameters for variuos models 
rq p b e carried out under a unified framework. Finally, the thesis directs attention 
to identify the topics for future research. 
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Chapter 1 

- 4 . 

Problem Sxa^ement 


1.1 Introduction 

In statistical signal processing we often assume the condition of stationarity or 
wide sense stationarity of the signal under investigation. Under the assumption of 
stationarity, people were interested in knowing the frequency domain characteristics 
of the signal. The periodogram along with different types of window are employed to 
find the power spectral density (PSD) of each signal [1], After that, the parametric 
modelling was introduced in determining the PSD. In parametric modelling, the 
signal is fitted into a known model of unknown parameters, and the analysis is 
done mainly on the behaviour of the chosen model rather than the signal. The 
main advantage of parametric modelling are flexibility of analysis and synthesis, 
and improved resolution in spectrum for short length of data. The improvement 
or resolution is due to the fact that we are providing more information about the 
signal. For example, in an autoregressive (AR) model for a finite length of data, one 
can compute only a finite number of auto-correlation function (ACF) from samples. 
But with the help of the AR coefficients, we can compute ACFs at any lag from the 
ACFs computed directly from the samples [1]. In this regard, it should be kept in 
mind that a bad choice of model can lead to a grossly erroneous result. 

It is to be stated, however, that a major weakness of the methods described above 
is the assumption of stationarity. In fact the real signals are often non-stationarity 
in nature, and we can name a few of those such as speech, ECG and EEG signals 
[2], The first attempt to tackle this problem is to introduce the notion of quasi- 
stationarity. Here the signal is assumed to be stationary only for a short length of 
time. This approach helps to justify the validity of the model for the signal. But 
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as the data length becomes shorter and shorter, the estimate of parameters looses 
its accuracy. So, one has to compromise between two contradictory requirements 
to get a workable model. In speech signal analysis, a segment of 25ms is considered 
to be good for most of the phonemes. But for the quickly varying phonemes such 
as plosives, the time duration is too large [3]. Hence, one has to compromise with 
the accuracy of parameter estimate. It seems that the stationary model cannot be 
expected to give good results when the signal under investigation is non-stationary 
in nature. 


1.2 Literature Survey 

1.2.1 Evolutionary Spectrum 

For a random process {X(t)} whose second order moments are finite and mean is 
zero, the covariance function R is defined by 

R(X-,r,t) = E{X(t)X(r)}, (1.1) 

where E is the expectation operator and (“) indicates complex conjugation. 

For a stationary process, the spectrum is the Fourier transform of the covariance 
function. The spectrum gives the information about the distribution of signal energy 
over the frequency axis. The idea of an evolutionary spectrum is just an extrapola- 
tion of the concept of spectrum from the stationary case to the non-stationary case. 
When the process is non-stationary, the covariance function not only depends upon 
the lag ( r — t ), but also on the instant r and t. So, it can be said that the evolution- 
ary spectrum will provide information about the distribution of signal energy over 
the frequency axis at each instant of time. That is in other words, the evolutionary 
spectr um is a measure of energy density over the time-frequency plane. 

In Reference [4], Loynes has given a set of properties, desirable or essential for 
spectrum of a process which is not restricted to be stationary. The list is as 
follows : 

A1 : The spectrum ’F is a real function of time and frequency, and is completely 
determined by the covariance function. 

A2 : T describes the distribution of energy over frequency. 
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A3 : *3/ transforms reasonably and preferably simply when the process {X(t)} 

is transformed linearly. In particular, a. knowledge of the spectr um of X 
determines the spectrum of the transformed process. 

A4 : The relation between the spectrum ’F and the covariance R is one to one. 

A5 : reduces to the ordinary spectrum, or some simple transform of it, if A is 

a stationary process. 

A6 : If the process is composed of a succession of stationary parts, say X\(t) if t < 

0 and X 2 (t) if t > 0, then the spectrum is also composed of the corresponding 
section of the stationary spectra. 

A7 : \E r is estimatable in principle, probably from one (infinite) length of record. 

A8 : d/ is the Fourier transform, or some related transform, of some apparently 
meaningful quantity. 

Any description of evolutionary spectrum may not satisfy all these properties but 
A2 and A3 seem to form the minimum set which is to be satisfied. Attempts have 
been made to define the evolutionary spectrum in various ways both parametrically 
and non-parametrically. Note that our discussion here will consider only the discrete 
time processes and their evolutionary spectra. 

1.2.2 Non— parametric Approach 

Similar to the case of stationary signals, for non-stationary signals also a number of 
non-parametric models were presented initially. The short time Fourier transform 
(STFT) [5] and Wavelet transform (WT) [6] are examples of this class. Both of these 
procedures can be looked as transformations from time domain to frequency domain 
of the windowed signal [7] [8]. One consequence is that the shape of the window 
has its effect on the energy distribution in time-frequency plane. The other is the 
fact that the time-frequency resolution is bounded by the uncertainty principle [9]. 

A familiar non-parametric time-frequency response is Wigner Distribution (WD) 
which has a large number of desirable properties [10], [11], [12], [13]. Its physi- 
cal interpretation as energy distribution at any point in the time-frequency plane, 
however questioned as the WD may assume negative value [10], [14], [15], [16], 
[17]. Similar objection is raised about its interpretation as energy density at a 
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point because tlie measure takes integration over the entire time axis [18]. More- 
over, as a quadratic time-frequency distribution it produces interference in case of 
multicomponent signals [10], [16], [17]. Though considerable amount of research 
has been done to reduce interference by choosing proper kernel, the attenuation in 
cross-terms comes only at the cost of broadening of WD signal terms [16], [17]. 

1.2.3 Parametric Approach 

In case of parametric modelling, at first the adaptive models are used for non— 
stationary signals. Here the first requirement is that the measure of non-stationarity, 
which may be the rate of variation of the auto-covariance function with time, must 
be less than the speed of adaptation. Under this assumption, the error in es tima tion 
of weight vector has two components. The first one called the weight vector noise 
is proportional to the speed of adaptation, and the second one called the weight 
vector lag is inversely proportional to the speed of adaptation [19]. Hence, a com- 
promise has to be made to keep the total error variance minimum. The conflicting 
requirements of the components of error limits the use of adaptive model in case of 
highly non-stationary signals. 

Next, in Prony’s model, the represented signals do not require to be stationary 
[1], This model, however, is very sensitive to noise. With the introduction of SVD 
decomposition [20] and auto-correlation like functions [21], [22], sufficient noise 
immunity is achieved in the method solved non-iteratively. Some iterative methods 
are also developed to estimate the parameters in the MLE sense [23]. The technique 
has got its application in modelling voiced speech signal [24], This model, however, 
is useful if the signal is a decaying signal. Otherwise, the model may require large 
order and the estimate of the model parameters will be difficult [1]. 

The approach of parametric modelling of non-stationary signals by the use of 
time— varying rational transfer function (RTF) model [25], [2] significantly does not 
have the short co ming s of the other methods, as mentioned earlier. Under certain 
simplified conditions, the parameters of the tune— varying RTF methods can be 
estimated, and the model has found application in case of natural signals like speech 
[26] and bat’s echo [27]. 
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' 1.3 Objective of Present Work 

A considerable amount of work has been done in parametric modelling of non— 
stationary signals. This thesis is devoted to the extension, consolidation and inte- 
gration of various approaches used in this area. The time-dependent RTF model is 
designed to estimate the evolutionary spectra of non— stationary signals [2], In the 
developed method, the main thrust is given to estimate the time— dependent AR 
coefficients of the fitted system. The issue of time— dependent gain is not properly 
addressed in this work. In nature, many signals like speech, have time-varying 
excitations in their production systems [28]. Hence, the use of time-varying gain 
in estimating the evolutionary spectrum becomes important. In Chapter 2 of this 
thesis, we have presented a non-iterative technique to estimate the time-varying 
gain parameter. 

It is to be noted that beside the time-dependent RTF model, several non- 
stationary parametric models are in use. These are Prony’s model, the amplitude 
modulation (AM) model, the frequency modulation (FM) model, and the amplitude 
frequency modutalion (AFM) model. So, it becomes an interesting issue to know 
whether these models are equivalent to the time-dependent RTF model, and if so, 
what are the conditions imposed on the general time-dependent RTF model for the 
model to be specific. Chapter 3 of this thesis is devoted to this specific issue of 
equivalence and unification. 

The electrocardiogram (ECG) is a real life non-st ationary signal. This signal 
has prominent pseudo-periodicity and burst like QRS complex makes it unfit for 
forced modelling into a stationary predictive model [29]. The various subsignals of 
ECG ( P, QRS and T ) represent different activities of heart [30]. So, from the 
modelling point of view, it is more meaningful to model these subsignals separately 
instead of the ECG signal as a whole. It is shown that the amplitude modulation 
model is a natural choice for this type of signal. Chapter 4 is devoted on the topic 
of analysis and synthesis of the ECG signal. 

To summarise the thesis, it can be stated the various aspects of parametric 
modelling of non-stationary signals are addressed here. The relationship exhibited 
between the time— dependent ARMA process and, each of AM, FM, AFM models 
and Prony’s model, provides important insight about non-stationary modelling. 
The thesis also identifies the problem for further research in this direction. 



Chapter 2 


Time-Dependent AR, Mode: with 
Time- Varying Gain 

2.1 Introduction 

Parametric modelling of non-stationary signal by employing the time— dependent 
rational transfer function (TRTF) model was proposed by Grenier [2]. The repre- 
sentation of time-varying coefficients by some basis functions gives this approach 
convenient form, both in terms of estimation and storage. Though the aspect of 
time variation of gain or changing variance of the noise input is considered at the 
time of formulation of the problem [2], its implementation together with a workable 
estimation technique is not found in literature. Even the simplest kind of non- 
stationary model where the AR coefficients are independent of time and only the 
gain parameter is varying with time, has not been studied to the best of our knowl- 
edge. In case of signals like speech, where the input energy of the production system 
is time varying [28], the time-varying gain becomes an important model parameter. 
It is therefore proposed here that the effect of time-varying gain in estimating the 
evolutionary spectra be studied with the help of both synthetic and natural signals. 

It is anticipated that the phonetic segments in speech waveform are characterised 
by their invariant spectral shapes. As a consequence, gain becomes the only pa- 
rameter that is time - 1 varying in nature in this case. The suitability of autoregressive 
model with time-varying gain only (ARTG) is studied for speech phonemes. 


6 



7 


2.2 Problem Formulation 

Let us consider a non— stationary sequence x(n)j n = 1,2, ,7V which is the output 

of a time-dependent autoregressive process with time-varying gain (TARTG) [2, 

31], 

x(n) + oi(n)x(n — 1) 4 h a P (n)x(n — p) = g(n)w(n) (2.1) 

for n — 1,2,... TV ; where g(n) is the time-varying gain, p is the order of the AR 
process, w(n) is the zero mean white noise input. 

Note that the gain g(n) does not have any effect on the locations of the poles. 
The pole locations are completely determined by the AR coefficients ai(n). The 
residues at the poles are decided by both a,(n) and g(n). Hence, to describe the 
system completely both 0,(12) and g(n) are to be computed. It should also be noted 
that the coefficient ao(n) of x(n) is set to unity in (2.1). This is conventionally 
done to obtain a non-trivial estimate for a = [ao(n) ai(n) • •• Op(n)Y while min- 
imising the noise variance of the AR process [32]. To get a non-iterative method 

for estimation of g(n ), let another sequence y(n) be considered such that, 

♦ 

/ N X ( n ) frs r\\ 

y{n) ~ g(n) (2 ' 2) 

where g(n) > 0 ,Vn, and n = 1 , 2 ,..., N. Combining Equations (2.1) and (2.2) we 
get, 

g{n)y(n) + a^n^n - l)y(n - 1) H + a p {n)g{n - p)y{n - p ) = g{n)w{n) (2.3) 

We now define a positive nonzero number e such that, 

e = sup 1 — ^ n . Vi = 1, ...,p. (2.4) 

n , 9(n) 

In case of real-life sampled signals, the e- factor can be made as small as desired by 
increasing the sampling frequency / s and as f$ — > oo, e — * 0. For such a signal, one 
can get from Equation (2.3), 

y(n) + ai(n)y(n - 1) H + Op{n)y{n - p) « w{n ) (2.5) 

which shows that the sequence y(n) is the output of a time-dependent autoregressive 
(TAR) process with constant gain. The model is of order p, and it is driven by the 
input, sequence u?(ti). ffow the time varying coefficients can be determined with 
the technique given in Appendix A [2]. 
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It remains to show how one gets y(n) from x(n), in other words, how the gain 
factor g{n) is estimated for each n. As it is already assumed that g[n) varies 
slowly with respect to time, one can estimate it as the local envelope or the short 
time energy of x(n), when the resulting sequence y(n) is of constant variance. By 
simulation it is verified that both the interpretations give similar results. Here we 
have chosen the first approach as it requires less computation. The estimate of g(n) 
will then be unique up to scaling. Hence, for a symmetric window ‘ 7 ’ of even length 
L, g(n) can be estimated as, 

1 L 

9 ( n ) = 7 £ M*)®(*» + * - §)1 (2.6) 

L »=i 

for n= (TV — £). 

The estimated envelope will be approximately unbiased only when the length 
of the window (L) is less than or equal to the order (p) of the AR model, and 
\Ya= 1 h'(i)y(n + i — |)| = 1 . The choice of L depends on the nature of the signal 
x(n) and that of the sequence g(n). For very small L, the averaging is not complete 
and the randomness of x(n) is reflected in g(n). On the other hand, if L is too 
large, it fails to follow the local variation of the envelope. As g(n) is only unique 

maxjx(i)| 

up to scaling, it may be scaled by a factor which makes g(n) suitable to be 

called the envelope of x(n). 

The non-stationary signal x(n) which is an output of a TARTG is modelled 
as another non-stationary signal y(Vi). amplitude modulated by g(n) as shown in 
following Figure 2 . 1 . The signal y(n) is the output of a TAR model with constant 
gain, and g(n) is an envelope function which slowly varies with time. 


TftRTG 



Figure 2.1: The block diagram of time-dependent AR model with time-varying 
gain (TARTG) 
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2.2.1 Evolutionary Spectrum 

The time— varying transfer function of the system, described by the difference equa- 
tion (2.1) can be represented as [33], 


H(n, u) = - 


9{n) 


ELi a k (n)exp(-j uk) 
and hence, the evolutionary spectrum of signal x(n) is given by 

P xx (n,oj) = \H(n,uj)\ 2 al 


(2.7) 


( 2 . 8 ) 


where cr is the variance of the white noise sequence w(n). 

It is known that the evolutionary spectrum of y(n) is given by 

Pyy(n,w) = ]Y*_ iak (n)expHu,k)f’ 
Hence, we find the relationship between the spectra to be 

P xx {n,u) = |s(n)| 2 P w (n,u;) 


(2.9) 


( 2 . 10 ) 


2.2.2 Spectral Distortion Measure 

A spectral distortion measure (SPDM) is a function of S and S which are the true 
and the estimated spectra respectively of the underlying process. The SPDM assigns 
a non-negative number d(S. S) to represent the distortion in using S to represent 
S. As the mean squared error does not appear to be subjectively meaningful in 
applications like speech coding, the new distortion measures have been introduced. 

A distortion measure proposed by Itakura and Saito, and termed as the error 
matching measure is defined as 

djs{S, S) = - ln(-^) — lj (2-11) 

o £> *1 

where S: and S are true and estimated spectrum of the estimated process and |[ • \\ p 
represents L p norm. 

Gray and Markel proposed a symmetrised Itakura— Saito distortion measure 
termed as the COSH measure of the form [34], 

dcosii{S, S) = ~[dis{S, S) + dis{S,S )] 


( 2 . 12 ) 
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For our computation, we have used the symmetrised Itakura— Salto distortion mea- 
sure or the COSH measure. To use it for the evolutionary spectra, the averaged 
spectral distance (ASPD) is introduced and it is defined as 

D(S, S) = - £ dcosHiSu St) (2.13) 

n i= 1 

where 5 = (Si, S 2 , ... , S n ) and 5 = (Si, S 2 , . . . , f? n ) are the true and the estimated 
evolutionary spectra respectively [34]. 


2.3 Simulation Results 

Example 2.1 

Consider an AR model of order p — 2, 

p 

x(n ) + ^ a.jx(n — i) = g(n)w(n ) 

i= 1 

where u)(n) = Gaussian white noise with zero mean and unit variance, and 

g(n) = l + 0.5sin( l ^ o ) 

The white noise w(n ) is generated with the help of a NAG subroutine. The model 
is started with zero initial condition, and the data of 5000 samples is recorded after 
4000 initial samples in order to get rid of the effect of transient. The estimate of 
g(n), g(n ) is done with the Hanning window of length L = 400 samples. The signal 
x(n) and its envelope g(n) are plotted in Figure 2.2. With the help of g(n), the 
normalised signal y(n ) is found, and the AR model is fitted into it. The statistics 
of coefficients for 100 simulations are given in Table 1, which roughly shows the 
validity of the model. The evolutionary spectra of the signal and its estimates with 
the TAR with constant gain and the TARTG are given in Figure 2.3. For the TAR 
with constant gain a sufficient number (up to 11) of Fourier basis functions are tried, 
but it is found that no order of basis can produce the time variation of evolutionary 
spectrum satisfactorily, the time variation being caused by variation in gain. Using 
the symmetrised Itakvra-Saito spectral distortion measure (SPDM), we found that 
the spectral distortion averaged over time (ASPD) of the estimated evolutionary 
spectra are 7.56X10 -1 and 2.09X10 -2 for the TAR with constant gain and TAR 
with time-varying gain respectively. 



11 


Table 2.1: Bias and variance in estimating the parameters by TARTG model in 
Example 2.1 


Serial Parameter Bias Variance 

i a i 

1 -0.5 3.0X10" 3 2.09X10" 4 

2 0.49 5.38-X10 -3 1.36X10~ 4 


D 

3 

+J 

■H 

C 

0 


sample number 

Figure 2.2: The original signal, and estimated envelope of the signal in Example 



2.1 
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Example 2.2 


In this example, an AR model of order p — 2 is considered. The TAR coefficients are 
chosen as ai(n) = —(0.5 — 0.002n), and < 22 (^) = 0.49. Note that the time— varying 
gain or standard deviation (SD) of white noise is linear in time. The function g(n) 
is given by 


g{n) = 1 + 


0.25n 

1000 


A set of 5000 samples are considered. The estimate of envelope g{n) is done with a 
Hanning window of length L = 400 samples. The signal and estimated envelope of 
the signal are shown in Figure 2.4. The evolutionary spectrum of the signal, and its 
estimates with the TAR with constant gain and the TAR with time-varying gain 
are given in Figure 2.5. For the TAR, 2 polynomial basis functions are used, and it 
is verified that more bases do not effect any appreciable change in the evolutionary 
spectrum. Using the SPDM, we found the ASPD of the estimated evolutionary 
spectra axe 1.877.XT0 -1 and 1.82X10 -2 for the TAR with constant gain and the 
TARTG respectively. 



sample number 


Figure 2.4: The original signal, and estimated envelope of the signal in Example 
2.2 
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Example 2.3 

The signal considered in this example is a sequence of 6000 s amp les of the vowel 
sound / o/ sampled at 16 KHz. The envelope of the signal is estimated by a Hanning 
window of length L = 400. The signal and envelope are plotted in Figure 2.6. 
Here the original evolutionary spectrum is approximated by the short-time AR 
spectrum using Hanning window of length 400 i.e. 25 ms. While using the Fourier 
basis, the TAR model order 10 and number of basis 7 are found to be s uffi cient for 
accurate modelling. The approximated original and estimated evolutionary spectra 
are shown in Figure 2.7. 
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Figure 2.6: The original signal, and estimated envelope of the signal /o / (vowel) in 
Example 2.3. 
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Example 2.4 

Th.e signal considered, in tliis example is a sequence of 4400 s am ples of a fricative 
sound /s/ sampled at 16 KHz. The envelope of the signal is estimated by a Hanning 
window of length L — 400. The signal and envelope are plotted in Figure 2.8. Once 
again, the original evolutionary spectrum is approximated by the short-time AR 
spectrum using Hanning window of length 400 i.e. 25 ms. The TAR, model order 10 
and the number of Fourier basis 7 are found to be sufficient for accurate modelling. 
The approximated and estimated evolutionary spectra are shown in Figure 2.9. 
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Figure 2.8: The original signal, and estimated envelope of the signal /s/ (fricative) 
in Example 2.4. 



18 



Figure 2.9: The evolutionary spectra of 
mated by (a) the QSAR model; (b) the 


(fricative) in Example 2.4, esti- 
^ del; (c) the TAR model. 
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2.3.1 Discussion 

In all the examples cited in this chapter, the estimated evolutionary spectrum con- 
sidering the time variation of gain is visually more akin to the original evolutionary 
spectrum. In Example 2.1, the ASPD of the evolutionary spectrum estimated by 
the TAR with constant gain is more than 30 times than the ASPD of evolution- 
ary spectrum estimated by the TARTG. In Example 2.2, this factor is 10 for the 
TAR with constant gain over the TARTG. This clearly shows the superiority of the 
proposed model in the synthetic examples. 

In case of natural signals like speech, the original evolutionary spectra are not 
known, and hence, quantitative tests are not possible in these cases. Nevertheless, 
we can see that the estimated evolutionary spectrum by the TARTG is more akin 
to the approximated evolutionary spectrum by short time AR spectrum than the 
similar estimate by the TAR with constant gain. This justifies the use of the TARTG 
in case of natural signals also. 

2.4 Modelling of Speech Phonemes by Autore- 
gressive Model with Time— varying Gain 

2.4.1 Introduction 

Speech is the oldest and most effective mode of communication in human society. 
In various speech processing problems e.g. coding, recognition, etc., parametric 
modelling of non— stationary speech signal is an important step. The conventional 
approach assumes stationarity for the signal over short time- frames, and estimates 
the set of parameters over each such frame [3]. The usefulness of this approach of 
modelling, which will be referred to as the quasi-stationary model in this study, be- 
comes limited because the evolutionary nature of signal is not revealed or exploited. 
This results in an inefficient modelling. The redundancy of information is termed as 
interframe redundancy and various coding techniques have been devised to exploit 
it [35]. 

On the other hand, the time-dependent rational transfer function model devel- 
oped in References [2, 26], presents a general approach for non-stationary signals. 
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The method, however, may become less attractive when a large number of basis 
functions are needed - to express the time-variations of the prediction coefficients, 
and the model order is sufficiently high. 

In this thesis, therefore, a middle approach of modelling based on phonetics is 
proposed for speech signals. Similar concepts are already in use for speech coding, 
where the need of variable frame length for the quasi- stationary model is advocated 
[36, 37]. The idea is to relate modelling to linguistic patterns like phonemes, and 
propose simple non-stationary model valid for each such speech segment. 

It is anticipated that phonetic segments in speech waveform are characterised 
by their invariant spectral shapes. As a consequence, the autoregressive model with 
time varying gain is proposed in this thesis. The usefulness of the proposed model 
is demonstrated by extensive simulation with a variety of phonemes. 

2.4.2 Model Building 

The autoregressive model with time-varying gain (ARTG) is only a special case of 
TARTG, where all the AR coefficients are time invariant. So, all the Equations 
(2.1) - (2.10) for TARTG are also valid for ARTG with corresponding change as 
mentioned above. 

Parameterisation of Time-varying gain 

When the time-varying gain g(n) is a slowly varying signal, it can be parameterised 
using a small set of appropriate basis functions. An expansion in terms of cosine 
functions is found suitable for g(n) and the cosine half series (Appendix B) is trun- 
cated where the first turning point of the estimated error variance occurs [38]. In 
general the change of error variance with the order of expansion is very slow in this 
study. Hence, an alternative criterion for selecting series-order may be set, which 
is to have the error variance less than (1/100) fraction of the signal power. If both 
of these strategies fail, the order is taken as a preset maximum value of 30 in this 
case. 

Order selection of Stochastic process 

The most commonly used order-selection criterion of an AR model is the Akaike 
Information Criterion (AIC). This criterion is known to be an inconsistent one, 
and it usually overdetermines the model order [39]. Hence, a consistent and more 



21 


generalised criterion, called the Efficient Detection Criterion(EDC), is chosen. For 
AR model, the EDC takes the form [40] : 

EDC (p) = AT log cr + C N p (2.14) 

such that 

hm ( C n /N ) = 0, 

A~*oo 

; lhn {C N /\og\ogN) = oo, 

where p, q, and N are as defined earlier. The penalty function is to be chosen 
according to the signal under study. In this work, Cv is set to be Cn = log N, 
for which EDC takes the form of Bayesian Information Criterion [41], Finally, the 
order p should be so chosen such that the EDC(p) is minimised. 

2.4.3 Results and Discussion 

For the purpose of analysis, a total of 65 English phonemes consisting of 8 cardi- 
nal vowels and 57 consonants, supplied by Central Institute of English and Foreign 
Languages, Hydrabad 500007, India, is taken into consideration. Each of the car- 
dinal vowels is uttered four times by the same speaker, out of which two sounds 
are of extended utterance. Each of the consonants is also uttered four times by 
three different male and one female speakers. The analog data supplied is low-pass 
filtered with cut off frequency 7.5 KHz and sampled at 16 KHz. As each phoneme 
gets mixed at both the ends by the following and preceding utterances, a segment 
of maximum length of uniform characteristic like the average rate of zero-crossing, 
is chosen by visual inspection of the waveform. To facilitate the comparison of the 
developed method with the quasi-stationary AR model (QSAR), the length of each 
segment is restricted to multiples of 25 ms i.e. 400 samples. 

With the data sets of different lengths in hand, the first step is to find the 
time-varying gain g(n). The rectangular window of length 10 ms consisting of 
160 samples, is found to be good for this purposed As the gain is unique up to 
scaling, it becomes meaningful to call it the envelope of the signal by introducing 
appropriate scaling factor.The envelope can be parameterised by the cosine half 
series. The orders of the cosine series axe calculated for each phoneme and the 
average is computed. The average number (= 18) of basis functions needed is 
highest for vowels, and it is almost double than the average number for consonants. 
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The average series-order over all speakers and all phonemes is 10. In Figure 2.10, 
it is shown how the choice of the number of basis functions affects the closeness of 
approximation of the envelope. The highest order of 30 provides the best result, 
whereas the order 10 produces too smooth an envelope. 
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Figure 2.10: The estimated envelope of the signal /o/, and its approximations for 
30,18 and 10 basis functions of cosine half series. 


Estimating first the gain g(n), the normalised sequence y(n) can be determined 
and then the corresponding AR parameters can be estimated. By using the EDC, 
the optimal order of the process is determined and the corresponding SNR is cal- 
culated. The SNR is defined as 


SNR = 10 log 


signal power 
prediction error power 


For the fricative /s/, the original signal and the signal regenerated by the ARTG 
are shown in Figure 2.11. For the ARTG, the average of optimal AR orders for all 
the phonemes is found to be 8, and the average SNR is calculated to be 26 dB. For 
the vowel sounds, the average SNR is 36 dB, and it is 24 dB for consonants. 

The merit of the ARTG model may be judged in relation to the other existing 
methods. For the quasi-stationary AR model (QSAR) the stationarity is assumed 
over a period of 25 ms i.e. for 400 samples [3]. For the non-stationary AR model 
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Figure 2.11: The signal /s/, (a) the original, and (b) the synthesised one by ARTG 
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(NSAR), the prediction coefficients are all functions of time [2, 26]. The order is 
chosen to be 10 in each case. In case of the NSAR model, 11 basis functions for 
each coefficients are employed, hence a total of (10 x 11 1 = 111) parameters are 

needed for each phoneme. For the QSAR model, the number of parameters will 
be (10 + 1 = 11) per block, where a block consists of 400 samples over a duration 
of 25 ms. For the ARTG model, the number of basis functions is chosen to be 15 
instead of the average basis-order of 10. This is necessary to get reasonably good 
approximation in case of vowels. Hence for the ARTG, a total of (10 + 15 + 1 = 26) 
parameters per phoneme is required. The achievable SNRs for all the models are 
tabulated in Table 2.2, whereas the number of parameters required per phoneme in 
each case is given in Table 2.3. Table 2.2 shows that the ARTG performs almo st 
same as the NSAR model (0.5 dB less) and little inferior (~ 2 dB) to the QSAR 
model. It can be shown from Table 2.3 that the QSAR model us uall y needs the 
maximum and the ARTG needs the minimum number of parameters to represent a 
phoneme. If the number of parameters per phoneme required by the ARTG is taken to 
be 1 p.u., then the QSAR and NSAR models need 4.95 p.u. and 4.27 p.u. respectively. 


2.5 Remarks 

The study reveals that if a non-stationary signal has a variation in energy, then 
the time-varying gain must be incorporated in modelling the signal. In this way, a 
faithful estimate of the evolutionary spectrum of the signal can be obtained. It is 
clear that time variation of the AR process cannot represent the effect of variation in 
gain. Though the proposed modification costs a little more computation, it makes 
the model more meaningful. In many processes of signal generation, the input signal 
energy is known to vary with time e.g. in speech signal production [28]. Recent 
study shows that the proper understanding of the evolutionary spectrum helps to 
enhance the performance of coders for speech and image signals to a great extent 
[42]. 

It is revealed in the study that the ARTG model is more suitable for speech 
signal than the QSAR model or TAR model with regard to the parsimony of param- 
eters.The SNR for the ARTG is similar to that of the TAR model, and marginally 
inferior to the SNR of the QSAR model. The evolutionary spectrum computed 
by the ARTG is similar to that by the QSAR model, suggesting that the ARTG 
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Table 2.2: SNR in modelling phonemes by QSAR, TAR and ARTG models 


Phoneme 



SNR 





fi ± a 





dB 


Name 

Number 

QSAR 

TAR 

ARTG 

Plosive 

13 

28.78 ± 9.22 

27.22 ±9.2 

26.17 ±8.95 

Nasal 

7 

30.36 ± 8.53 

28.35 ± 8.53 

28.30 ±8.83 

Trill 

2 

24.51 ± 7.3 

23.35 ± 7.88 

21.49 ± 5.4 

Tap or Flap 

2 

27.97 ± 8.49 

27.01 ± 8.76 

26.73 ± 8.62 

Fricative 

22 

23.28 ± 8.49 

21.96 ± 8.37 

21.22 ± 8.17 

Lateral Fricative 

2 

23.26 ± 7.19 

22.3 ±7.48 

20.72 ±6.79 

Approximant 

5 

29.13 ± 9.52 

27.22 ± 9.45 

27.52 ±9.73 

Lateral 

4 

30.15 ± 8.89 

28.31 ± 8.93 

27.43 ± 8.91 

Approximant 





Vowel 

8 

37.8 ±5.06 

35.12 ±4.88 

35.84 ±5.44 

Average 


27.98 

26.32 

25.8 



26 


Table 2.3: Number of parameters required by QSAR, TAR and ARTG models for 
phonemes 

Phoneme Number of Parameters per Phoneme 


Name 

Duration 

H ± <7 

ms 

QSAR 

TAR 

ARTG 

Plosive 

246.9 ±147.43 ■ 

108.64 

111 ' 

26 

Nasal 

243.83 ± 88.73 

107.28 

111 

26 

Trill 

268.75 ± 124.84 

118.25 

111 

26 

Tap or Flap 

184.38 ± 128.05 

81.13 

111 

26 

Fricative 

264.77 ± 150.64 

116.5 

111 

26 

Lateral Fricative 

218.75 ± 93.33 

96.25 

111 

26 

Approximant 

277.5 ± 143.81 

122.1 

111 

26 

Lateral 

Approximant 

256.25 ± 142.11 

112.75 

111 

26 

Vowel 

564.84 ±246.6 

248.53 

111 

26 

Average 

292.56 

128.72 

, 111 

26 
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is the best choice to model phonemes. The entire study suggests that within a 
phoneme the change in spectral wave shape is small whereas the short-time energy 
can have a large variation. Figures 2.7 and 2.9 also support the above-mentioned 
fact. Moreover, the time-varying gain or the envelop estimated as a part of the 
model identification, can be of use in the recognition of consonants [43, 44]. There- 
fore, the ARTG seems to be quite promising in modelling phonemes, which can in 
turn be used directly in speech synthesis. Moreover, the ARTG model together with 
a phoneme boundary detection algorithm can be used for coding and recognition of 
the continuous speech. 



Chapter 3 


Modeling Xon— Stationary Signals 
: A Unified Approach 


3.1 Introduction 

Under the general framework of time-dependent ARMA (TARMA) model, any sig- 
nal can be represented as the output of a system with time-varying poles. It is 
revealed in this thesis that an interesting classification of the non-stationary sig- 
nals is possible through functional specifications of the movements of time-varying 
poles. In this regard, this study presents a unified approach to encompass various 
techniques of parametric modelling of non-stationary signals viz. Prony’s model, 
the amplitude modulation (AM), frequency modulation (FM) and, amplitude and 
frequency modulation (AFM) models [21, 45, 46, 47, 48]. 

In my study, it is found that the damped sinusoid, AM, FM and AFM sinu- 
soids can be generated from the TARMA model by using proper restrictions on 
the model parameters. This finding provides a general procedure for estimating 
the parameters of various non-stationary models using a set of basis functions. In 
this way, the estimation of parameters can be obtained by solving a set of linear 
equations as shown in Appendix A [2], It is demonstrated that the unified approach 
of signal modelling and parameter estimation results in greater modelling flexibility 
(choice of more suitable model) and/or -accuracy. The better accuracy is achieved 
because the method presented here transforms the nonparametric/nonlinear esti- 
mation problems into the respective parametric/linear estimation problems. 

This thesis reveals that the impulse driven time-dependent all pole models viz. 
the AM, FM, AFM and dam ped sinusoidal signal models can be visualised as some 
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Figure 3.1: Alternative non-stationary models 

special cases of the time-dependent AR (TAR) model. This fact is explained in the 
following chapter and the interrelation among various models is shown in Figure 
3.1. 

The trajectory of pole of the impulse driven system keeps information of the 
modulating signal in case of the modulating signal model. For the damped sinusoidal 
model, the pole trajectory provides information about the rate of decay of the signal. 
Note that the direct estimation of poles gives rise to a non-linear estimation problem 
even in case of stationary signals. The situation becomes much more complicated 
with non-stationary signals. 

It is to be pointed out that although the TAR coefficients contain the same 
information about the modelled signal, the coefficients do not provide any charac- 
terisation for the signals. But the TAR coefficients can be estimated by solving a 
set of linear equations, when the coefficients are expressed in terms of some suitable 
basis functions. This motivates us to estimate the TAR coefficients and then, to 
find the pole trajectories by solving the roots of the polynomial. In this way, we 
get a meaningful characterisation of the signals through a relatively straightforward 
estimation technique. 
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3.2 Model Building 


3.2.1 Problem Formulation 

The input /output relation of the process is given by [31, 49], 


x(n) + ai(ri)x(n -!) + ••• + a r (n)x(n - r) = b 0 (n)w{n) + • • • + b q {n)w{n - q) (3.1) 


where a k (n) : k = l,...,r are the AR parameters, bi(n) : l = 0,1,..., g are the 
MA parameters, w(n) is the zero mean white noise input, and x(n) is the output 
sequence of the process. The sequence x(n) is assumed to be non-stationary. 

The generalised transfer function H ( n , z) corresponding to the diff erence equa- 
tion (3.1) can be expressed as [49], 


H(n,z) 


Ef=o h(n)z~ l 


(3.2) 


where it is assumed that for all n, 


H(n,z) ~ H(n — k,z ), for k = 1,2, . . . ,r (3.3) 

The condition imposed by (3.3) is generally valid for slowly varying system, and 
it is usually less restrictive than the assumption of stationarity for x(n) over a data- 
window. Note that, the model order is usually much small compared to the window 
size in case of the quasi-stationary AR model or the short time Fourier transform. 

The time-varying rational transfer function in (3.2) can be rewritten as [49], 


rr, . 9(1)11?.! (l-*(n)z- 1 ) 

H{n ' z) n 


(3.4) 


or, 


H(n,z) 


y- cjt(n) 
^1 -p k {n)z~ x 


(3.5) 


where g{n) is the time— varying gain of the generalised transfer function, Pk{ n ) '■ k = 
1, . . . , r are the time-varying poles, zi{n) : Z = 0, 1, .... g are the time-varying zeros, 
c fc (n) are the time-varying residues at the poles Pk( n ), and r > q for realisability. 

While writing (3.5), it is assumed that the rational transfer function has no 
repeated pole and the order of zero is less than the order of pole. Under these 
conditions, a rational transfer function can be written as a weighted sum of single- 
pole systems, where weights may vary with time. 
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It will, therefore, be instructive to study the first order time— varying AR mo del 
and understand the multiple-pole system as a combination of single-pole systems. 
Let a first order AR process be given by 


x ( n -) _ [ ax ( n ~ 1) + w i n ) if n < no 
w ad(n)x(n — 1) + w(n) if n > no 


| 3 . 6 ) 


where w(n) is the white noise input with zero mean and var ian ce 

Note that initially the coefficient ‘a’ is set at a = (1.0 - e) exp {j6) with 
e and the variance of the white noise is fixed at cr w . The values of e and o^ a 
made sufficiently small so that x(n) becomes a sinusoidal signal before no. This 
feature is explained in Appendix C [1]. At the time n = n o with no large enoug| 
to assure complete decay of transients, the time-varying factor d(n) is introduced 
and its effect on x(n) is investigated. 

In (3.6) as ja| — > 1 and cr^ — >■ 0, the time series x(n) grows to a sinusoid with 
amplitude much bigger than cr w [1], So, for n > n 0 , the contribution of w{n) in ±{»^| 
is negligible compared to the contribution of x(n - 1). Then the skeleton structure 
[50, Chapter 4] of (3.6) becomes 

x(n) = ad(n)x(n — 1), n > no (3.7| 


At n = n 0 + m and with repeated substitution leads to 

no-bm 

x(no + m) = a m d(l)x(n o) 

J==tio+l 

= a m D(m)x(no) 


(3.8) 

(3.9) 


where D(m) = n^i d(l) 


3.2.2 Classification of models 

It is shown in this section that the use of some mathematical functions for d(n) may 
lead to a few well known models for non-stationary signals. These models are, 

• Damped sinusoidal model, 

• Amplitude modulated sinusoidal model, 

• Frequency modulated sinusoidal model, 

• Amplitude and Frequency modulated sinusoidal model. 

In the following subsections, these models are investigated in detail. 
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Damped Sinusoidal Model 

In this case, the time-varying coefficient is set to be 

d(n) = 6 > 1, n > no (3.10) 

Substituting d(n) in (3.9) and after simplification we get 

D(m) = b~ 2 ; 

and 

x(no + m) ~ x(nQ) exp ( j6m ) (3.11) 

where a « 1.0 exp (j0). 

So, for d(ri) as defined in (3.10), the process decays rapidly as indicated by 
(3.11), and after sometime, the output becomes same as the input noise sequence. 
Note here that as a special case, if we set 


d(n) = b l , b > 1 & n >._ n 0 


(3.12) 


from (3.9) we get 

D(m) = b~ m 

and 

x(n 0 + m) ~ 6 _m s(no) exp (jdm) (3.13) 

where a ~ 1.0 exp (J9). The system loses its energy exponentially and settles to a 
new AR model with pole at l ab\ 

The difference between d(n) in (3.10) and (3.12) is that the damping factor given 
by In |d(n)| decreases linearly in the first case and it is constant in the second case. 
This not only results in rapid decay of the signal in the first case, the rate of decay 
also increases with time. 


Amplitude Modulated Sinusoidal Model 


In this case, the time-varying coefficient is substituted as 

d(n) = exp ( bslna(n — no)), n > no (3-14) 

which after substitution in (3.9) and with subsequent simplification leads to 

no-bm 

D{m) — IQ d{l) 

]=n o+l 


exp 


' ' — — (cos am — 1) 


> 

} 




(3.15) 
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The approximation of the above equation holds true for small a (Appendix D). 
Rewritten, 


x{uq + m) ~ 1.0 exp {jQm) exp (— (1 - cosam)^. x(n 0 ) 

V a J 


where a ~ 1.0 exp {jd) is substituted. 

Taking real part of the above equation, one gets 


real[x(no + m)] ~ A exp 


— (1 — cos am) cos (8m + 6) 
^a j ' 


(3.16) 


(3.17) 


where x(no) = A exp {jo). Note that (3.17) gives the output as an amplitude 
> modulated signal with modulating factor given by exp (—£ cos am). 

For any positive discrete-time signal ■u(n), we can modulate with the factor 


d(n ) = 


u (n) 
u(n — 1) 


(3.18) 


and, we get the desired output in discrete time domain as given by 


x(n 0 + m) = u(m) exp [j8m)x(n 0 )/u(n 0 ) 


(3.19) 


The time-varying part of the pole d(n) is real in the case of amplitude modula- 
tion. This results in time variation of the magnitude of the pole, keeping its angle 
constant. In other words, for a « 1.0 exp {jd), the pole will have radial movement 
alternating through the point a ~ 1.0 exp {jd), in this case. 


Frequency Modulated Sinusoidal Model 
In this case we assign the factor d(n) to be 

d(n) = exp (jbsina(n — n 0 )), n > no (3.20) 

which after substitution in (3.9) and with subsequent simplification leads to 

no+tn 

D(m) = J! d(l) 

l=nc+l 

~ exp ^ -j— (cos am — 1) (3-21) 

The approximation of the above equation holds true for small a (Appendix D). 
From (3.21) we obtain 

x(n 0 + m) ~ 1.0 exp (jdm) exp ^ (1 - cos am)\ x(no) 


(3.22) 
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where a « 1.0 exp (j9) is substituted as before. 

The real signal is given by 

real[x(no + m)}~ A cos [dm - - cos am + 6+ — \ f'3.23'1 

^ a a) ■ K J 

where ®(no) = A exp (j 6). It is clear that (3.23) gives the singletone frequency 
modulated output of the process. 

Therefore, in order to modulate with a discrete-time signal u(n), one has to use 


d(n) = exp ( jKu(n )) (3.24) 

for FM, which leads to the FM signal 

f 7io+m \ 

x(n 0 + m) = Aexp ' j(9m + K ^ « to + « , 

V i=no+l ) 


and 

d(n) = exp (jKu(n)) (3.25) 

for PM, where u(n) is the first order backward difference of u(n). The corresponding 
PM signal is 

x{n 0 + m) = .A exp ( j(9m + Ku(n) +(f> — Ku(n o)). 

The time-varying part of the pole d(n) is an imaginary quantity in the case of 
frequency modulation. This means that d{n) will cause only the angular movement 
of the poles. In other words, for a « 1.0 exp (J9), the poles will have circumferential 
movement alternating around the point a ~ 1.0 exp (J9). 

Amplitude and Frequency Modulated Sinusoidal Model 

In this case, we apply both the earlier mentioned schemes combinedly. The proposed 
d(n) will be of the form, 

d(n) = “ l(n) „<*p(jU 2 (»)) (3-26) 

. ' ui(n - 1) 

where the positive discrete-time signal Ui(n) is amplitude modulating signal and 
U 2 (n) is the frequency modulating signal. 
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3.2.3 Estimation of Parameters 

As a complex sinusoid can be represented by a first order AH, model, its real part 
will require a second order AR model. Let us have the characteristic equation 
corresponding to a first order time-dependent complex AR model (3.6) for n > 0; 

1 - a(n)z -1 = 0 (3.27) 

where a(n ) = ad(n) and z -1 is the delay operator [25, 51] . 

The characteristic equation of the corresponding real second order AR model is 
then 

(l — a(n)z -1 j^l — a*(n)z _1 ) = 0 (3.28) 

or, 1 — (a(n) +a*(n)) z _1 + a(n)z -1 a*(n)z -1 = 0 

or, 1 — (a(n) + a’(n)) z -1 + a(n)a*(n — l)z -2 = 0 

or, 1 — (a(n) + a*(n)) z _1 + |a(n)j 2 z~ 2 ~ 0 

where a*(n) is the complex conjugate of a(n ) and the system is assumed to be 
slowly varying in time as in (3.3). Let the second order time-dependent AR model 
to represent a real sinusoid be, 

x(n ) + ai(n)x(n — 1) + a 2 (n)s(n — 2) = w(n ) (3.29) 

where w(n) is the prediction error which takes into account all the sources of errors 
viz. the noise in the data, error in approximating the time-dependent prediction 
coefficients, etc. The sequence w(n) is assumed to be of zero mean and variance cr^,. 
Comparing (3.28) with (3.29) we get 

ai(n) = — (a(n) + a,*(n)) (3.30) 

= — 2|a(n)j cos (arg a(n)) 

= — 2|d(n)| cos (9 + arg d(n)) 

02 (n) ~ | a(n)] 2 (3.31) 

~ \d(n)\ 2 

where a « 1.0 exp (jd) and d(n) = |d(n)| exp (jarg d( n)). 

From (3.30) and (3.31) we get, 


\d(n)\ ~ \fa 2 (n) 


(3.32) 
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arg d(n) ~ arccos ^ - Q (3 33) 

\2sJaz{n)j 

From (3.18) the estimate of the modulating signal for amplitude modulation can be 
expressed as 

no+m 

wi(no + m)= I] \d{l)\u x (no) (3.34) 

i=no+l 

and from (3.24) the estimate of the modulating signal for frequency modulation can 
be given as 

&2 (n) = arg d(n) (3.35) 

Both these estimates give only the scaled and shifted version of the modulating 
signals. Note that \d(n)\ represents the radial movement of the system pole and 
the quantity keeps information about the amplitude modulating signal. Whereas, 
arg d(n) represents the circumferential movement of the pole and it keeps informa- 
tion of the frequency modulating signals. 

Since we can only get some scaled and shifted estimates of the modulating 
signals, as given in (3.34) and (3.35), we have to find the proper values for scaling 
and shifting. If u(n ) be the modulating signal, and u(n) be its estimate, for n = 
1,2 , . . . , IV, we have to find suitable values for the shift ‘p’ and scale l q' such that 
the expression En=i [u(n) — (p + qu(n))] 2 is minimised. The expression for p and q 
are of the form 


N (ZLl“( n M n )) - (En=l «(«)) (En=l U ( n )) 

&(»)*) “( Eir^n)) 2 


(3.36) 


E«=i u(n)-q E2-i u{n) 


r N 

The error variance of estimation is defined as 

N 


(3.37) 


- 1 £ [u{n) - {p + gfi(n))] 2 

-tV ““ 1 n= 1 


Essentially we estimate the modulating signal via the estimation of a.\ (n.) and 
a 2 (n). This is done by introducing a suitable set of basis functions for the coef- 
ficients, thereby, transforming a non-lineai estimation problem into a linear one. 
The details are given in Appendix A. 
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As the estimate from the parametric method requires scaling and shifting, these 
are also utilised for the estimates from the non-parametric methods which are de- 
scribed in the following paragraphs. In this way, we get a fair comparison among 
their performances of various methods for estimating the modulating signals. 


3.2.4 Non— linear Energy Operator 

For both continuous and discrete time signals, Kaiser has defined a non-linear 
energy tracking operator '3/ [52, 53]. For the discrete time case, the energy operator 
for x(n ) is defined as, 



\&[#(n)] = x 2 (n) — x(n — l)a;(n + 1) 

(3.38) 

For the signal 


x(n) = r(n) cos (<£(n)) 

(3.39) 

we have 


'^(n)] « (r(n)Q(n)) 2 
^/'^(n)] w |r(n)Q(n)| 

A 

(3.40) 


where f2(n) is defined as the instantaneous frequency i.e., Q(n) = The differ- 
entiation operator is approximated with a proper (backward, forward or symmetric) 
difference operator in the discrete time case. 


When one of the variables r(n) and O(n) is constant, we can get the other 
variable with a scaling of ^?[cc(n)]. So, the energy operator can estimate the 
modulating signal, or more precisely its scaled version, when either just AM or FM 
is present [54], 

When both AM and FM are present simultaneously, three algorithms are de- 
scribed to est ima te r(n) and f2(n) separately in [55]. The best among the three 
algorithms according to performance is called the discrete energy separation algo- 
rithm 1 (DESA-1) [55]. The DESA-1 is defined as follows 


y{n) = x(n) - x(n - 1) 


arccos 


f.Sb/Ml + »b(n + i)P 

V 4tf[i(n)] } 


fi(n) 


(3.41) 


^[x(n)] 



r(n) 


(3.42) 
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The frequency estimate works so long as 0 < O(n) < 7 r. 

3.2.5 Hilbert Transform Separation Algorithm 

The Hilbert transform of any AFM signal x(n) = r(n) cos ( 4>(n )) is a signal x(n) 
with the Fourier transform X(w) such that 

X(w) = -jsgn(u;)X(w) (3.43) 

where X(lj) is the Fourier transform of x(n). 

Given the analytic signal 

x(n) + jx(n) = s(n) exp (j9(n)) (3.44) 

its modulus s(n) and phase 9{n) can serve as (generally approximate) estimates of 
the amplitude and instantaneous frequency of x(n). Thus the Hilbert transform 
separation algorithm (HSTA) is given by [56] 

s(n) = \[x 2 (n) + x 2 (n) « |r(n)[ (3.45) 

9{n) = arctan « <f>{n) (3.46) 

From the instantaneous phase 9{n) we can get the estimate of instantaneous fre- 
quency by phase unwarping followed up with differentiation or (forward, backward 
or symmetric) difference in discrete domain. 

In practice, Hilbert transform can be implemented by using an FIR or a DFT 
approximation to the ideal IIR filter [57]. This approximation, however, introduces 
additional error in the Hilbert transform. The error increases with the decrease 
in the FIR filter length [56]. For this reason, we have implemented the Hilbert 
transform with a DFT over the full data length. 

3.3 Simulation 

3.3.1 Synthesis 

For the purpose of simulation we have taken n-o = 5000 so that the signal becomes 
steady within that t im e. The origin of the tune axis is shifted to 5000 for the 
purpose of displaying the output. In all the cases, |a| = 0.99999 and — 1.1 are 
used for simulation. 
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For generating the damped sinusoid with d(m) according to (3.10), we have set 

d{m) = (1.005) -m , 9 = -■ 

3 

and the real part of the signal as given by (3.11) is plotted in Figure 3.2. In the 
special case, we set d(m) according to (3.12), 

d(m) = 0.95 (a constant), and 9 = — , 

3 

and the output sinusoid settles to an AR process with a pole at 0.95 Z- and the real 
part of output as given by (3.13) is plotted in Figure 3.3. 

To generate amplitude modulated signal with <i(n.) as (3.14), we have chosen 

d(m) = exp (o.5 sin ) , B = | 

and the output real[x(no + m)] as given by (3.17) is plotted in Figure 3.4. As 
the d{m) is the exponential of a singletone, the process is named as exponential 
modulation. 

The singletone frequency modulated output is generated with d{m) as (3.20), 
d(m) = exp 0.4 sin B = | 

and the output real[x(n Q + m)] as given by (3.23) is drawn in Figure 3.5. 

3.3.2 Analysis 

Damped Sinusoidal Model 

We have generated two damped sinusoids according to (3.9) with : 

Case I : d(m) = 1.0 — 2.0 x 10 _5 m, 9 = 7r/6. 

Case II: d(m) = (1.00001)-*, 9 = ir/6. 

for 77 i = 1, ... , 512. The real part of the two sinusoids are shown in Figure 3.6 and 
3.7. In both the cases, we have fitted exponentially damped sinusoids according to 
Prony’s model and the time-dependent AR (TAR) model. For Prony’s model, the 
optimal model orders according to the EDC (2.14) are found to be nine and two in 
Case I & II respectively. For TAR model, two monomial bases are used in Case I 
and two Chebyshev bases are employed in Case H. The signal to prediction error 
ratios are computed for all combinations and shown in Table 3.1 . 
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Figure 3.2: The real part of generated damped sinusoid (rapid decay), Re[x(m)] 
= n£i d(l) cos (7rm/3), m = 1, . . . , 200; and d(l) = (1.005) -1 . 



Figure 3.3: The real part of generated damped sinusoid (transient), Re[x(m)] 
= n£i d(l) cos (■ ktu/S ), m = 1, ... , 200; and d(l) = 0.95. 
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Figure 3.4: The real part of generated amplitude modulated sinusoid (exponen- 
tial modulation) Re[x(m)] = Iliii d(Z) cos (7rr7T,/3), m = 1,...,200; for d(l) = 
exp (0.5 sin(^Z)). 
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Figure 3.5: The real part of generated frequency modulated sinusoid (single- 
tone modulation) i?e[x(m)] = Iliii ^(0 cos ( 7rm /^)’ m ~ 1, ...,200, and d(l) 
exp 0*0.4 sin (j^Z)). 
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Figure 3.6: The real part of damped sinusoid Re[x(m )] = ITHo d(Z) cos (mm/6), 
m = 1, . . . , 512; and d(l) = (1.0 - 2.0 x 10" 5 Z). 
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The real part of damped sinusoid iie[x(m)] — IIjILo cos (mm/6), 
512; for d(l) = (1.00001)-'. 
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Table 3.1: Signal to prediction error ratios in mode lling damped sinusoids. 

Case Prony’s Model TAR 
(dB) (dB) 

I 54.5 66 

H 58 77 


Amplitude Modulated Sinusoidal Model 

Consider the discrete time AM signal with singletone modulation, as given by 
y(n) = (1 + k cos u a n ) cos (9ri), n = 1, . . . , 512; 

where k = 0.8, oj a = 7r/128 and 6 = tt/ 6. The modulated signal is shown in Figure 

3.8. The modulating signal is estimated as given in Section 3.2.3. While using the 
Fourier and Chebyshev bases, 33 and 28 basis functions respectively are found to 
give best results. The estimated modulating signal using the time-dependent AR 
model (TAR) with Chebyshev bases and the error in estimation are shown in Figure 

3.9. 

The performance is compared to the results obtained using non-paramecric ap- 
proaches like the non-linear energy operator and Hilbert transform. The DESA-1, 
which provides the best performence among the three algorithms based on non- 
linear energy operators, is used. For the Hilbert transform implementation, a DFT 
over the full data length is employed for most accurate result [56] .The error variance 
in estimating the amplitude modulating signal by these methods are given m Table 
3.2. The loci of the poles are shown in Figure 3.10, which nicely shows the radial 

movement of the same on the unit circle. 

The estimation of AM signal by the proposed method is sensitive to noise. The 
Figure 3.11 shows the estimates of a 2 (n), and the corresponding estimations of 
the amplitude modulating signals are shown in Figure 3.12, when additive white 
noise is added to the signal. In order to introduce robustness in the estimation 
procedure, the system is over-determined using L = 7 in (A.12) (Appendix A). In 
this way, much noise immunity is achieved by considering the innovation at any 
instant orthogonal to the last L samples or more precisely [Y(n,l),...,r(n,p)], 
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Figure 3.8: The AM signal y(n) = (1 + fccosw a n) cos (0n), n = 1, , 512; and 

k = 0.8, u a = 7t/128, 6 = 7r/6. 


sample number 

Figure 3.9: The estimated amplitude modulating signal (solid line) and the estima- 
tion error (xlO) (dotted line). 
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Figure 3.10: The loci of the poles of the AM signal. 


Table 3.2: Error variance in estimating the amplitude modulating signal 
case. 

Non-Parametric Approach TAR Process 

implemented with 


Non-linear Energy Hilbert 
Operator Transform 


Fourier Chebyshev 

Basis Basis 


1.94 x 10~ 4 


3.7 xlO' 3 1.39 xlO- 4 4.65 X 10" 5 


in no noise 
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Figure 3.11: The estimates of the time-varying AE coefficient a 2 (n) for the AM 
signal at different SNR levels 



Figure 3.12: The estimates of the modulating signal for the AM signal at different 
SNR levels 
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when L>p, the AR model order. The additive white Gaussian noise is used for 
simulation, and 100 realisations are implemented with different noise sequences at 
each SNR level. The performance of the parametric method in comparison with 
the performance of the non-parametric methods are shown in Table 3.3. Simulation 
results for SNR below 40 dB are not shown because in such cases the energy of 
modulating signal occasionally becomes negative, leading to a meaningless result. 
This phenomenon is observed while using each of the methods. 

Frequency Modulated Sinusoidal Model 

Consider the discrete FM signal with singletone modulation, 

y(n) =cos{9n + uj m £? =1 cos (w/i)), n = 1,...,512; 

where u m = 0.26, uif = 7r/128 and 6 = 7r/6. The signal is shown in Figure 3.13. 
While using the Fourier and Chebyshev bases, 15 basis functions give best results 
in both the cases. The estimated FM signal in no noise case using Chebyshev basis 
and the error in estimation are shown in Figure 3.14. 

The method for estimating the frequency modulating signal is found to be robust 
when Gaussian white noise is added to the modulated signal. So, there is no need 
to make the linear system over-determined in this case. The locus of the poles is 
shown in Figure 3,15. 

The performance of the proposed method in comparison with that of the non- 
parametric methods are shown in Table 3.4. The results tabulated are obtained by 
100 realisations with independent noise sequences at each SNR level. The (*) symbol 
in Table 3.4 signifies the fact that the energy of modulating signal occasionally 
becomes negative in some realisations. 



Table 3.3: Error variance in estimating the amplitude modulating signal with 
overdetermined system. 
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Figure 3.13: The FM signal y (n) = cos(0n + u m E”=i cos (w,*)), n = 1, . . . , 512; and 
= = : 0.20, Cdf = 7t/128, 6 = 7r/6. 
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Figure 3.14: The estimated frequency modulating signed (solid line) and the esti- 
mation error (xlO) (dotted line). 
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Figure 3.15: The loci of the poles of the FM signal. 



Table 3.4; Error variance in estimating the frequency modulating signal. 


lO 


w £ 

CO £ 
CD 

cj -a 

o o 

Ch 44 

a 5 

B5 2 

H cd 





OO 

1 

oo 

x*- 

cn 

* p-H 



O 

rH 

o 

tH 

1 

O 

rH 

3 



X 

X 

1 X 

CQ 


t 

Tt« 

CM 


> 

<D 

fH 

CO 

IT 

4j 

o 

rH 

X 

q 

rH 

-H 

icq 

4 

41 

CO 

cm. 

+| 

£ 
0) 
rH 
— H 

3 

N* 

oo 

x- 

i 

O 

rH 

x-~ 

1 

O 

rH 

CO 

1 

O 

rH 

o 



X 

X 

X 




00 

LO 

q 

00 

rH 




00 

cS 

rH 


CQ b 

H 4J 

0) n 


c£ 


X 

CO 

00 


i 

O 

rH 

X 

CM 

Q0 

-H 


x 

CO 

q 

cd 


QO 

o 


00 

lO 


4) 


I 

o 


IN 

LO 

00 


oo 

i 

o 

1—1 

X 

CO 

CO 

ci 

41 


X 

b 


co ¥ 

o o 


o 

cd 

0 

Jh 

Oh 

CD 

C 

CJ 

L~l 

4-0 

CD 

a 

cd 

& 

Oh 

1 

d 

o 

2 


«H 

cS 

d ^ 
$ b 

& 4i 

44 53. 

fH 

CD 

X 


s 


J-t 

o 

44) 

cd 

fH 

<D 

a 

O 

hO ' 

u 

d 

W 

3 

CD 

d 


b 

41 

=1 


<0 

I 

O 


CO 

CO 


X 

q 

rH 

HH 

CO 

I 

o 


CO 

CO 


X 

10 

tq 

41 

co 

o 


CO 

CO 


X 


oo 


-H 


x 

oo 

CO 


d 

5§ 



t*. 

i 

<0 

LO 


O 

o 

1 

o 


rH 

rH 

rH 


X 

X 

X 

CO 

CM 

00 

N- 

o 

q 

q 

oq 

rH 

H 

rH 

rH 

X 

41 

41 

41 

lO 

CO 

L0 

H 

q 

b 

o 

1 

o 

CM* 

rH 

rH 

rH 


X 

X 

X 


CO 

00 

1> 


cq 

q 

o 


b 

cm 

cm 


CO 


CQ 


o 

CO 


o 

lO 


o 


CO 

t 

LO 


© 

rH 

b 

H 

I 

o 

rH 

X 

X 

X 

rH 

CO 

oo 

© 

CM 

cq 

cm* 


00 

41 

41 

4 

CO 

1 

LO 

H 

«— s 

b 

1 

o 

rH 

rH 

rH 

X 

X 

X 

lO 


rH 

LO 

o 

.cq 

cd 

d 

4 

X- 

LO 

LO 

b 

J 

<0 

b 

rH 

rH 

tH 

X 

X 

X 

L0 

L— 

CO 

N-’ 

rH 

rH 

+1 

CO 

41 

LO 

cd 

4 

1 

o 

1 

o 

H 

1 

rH 

rH 

o 

X 

X 

X 

© 

b- 

CO 

q 

q 

rH 

cd 

rH 


LO 

LO 

1 

LO 

b 

rH 

o 

rH 

1 

o 

rH 

X 

X 

X 

LO 

00 

H 

q 


100 

id 

-H 

41 

41 

CO 

1 

CO 

co 

O 

rH 

b 

rH 

b 

rH 

X 

X 

X 

C0 

q 

i>- 

cq 

CM 

rH 

cd 


* 


* 

LO 

lO 

lO 

o 

rH 

o 

rH 

o 

rH 

X 

X 

X 

© 

IH 

© 

o 

JH 

N- 

Q0 

© 

id 

4 

41 

41 

C0 

1 

CO 

CO 

O 

rH 

1 

O 

rH 

b 

rH 

X 

X 

X 

t- 

q 

,TJJ 

q 

rH 

4 

id 

o 

o 

lO 

CO 

CM 

rH 


central library 

I L T., KANPUR 

ta W* A. JJJ M5 


52 


Amplitude and Frequency Modulated Sinusoidal Model 
Consider the discrete AFM signal, 

y{n) = (1 + k cos u a n) cos ( 9n + cos (w/i)), n = 1, . . . , 512; 

where, k 0.5, bj m 0.50, ui a — ujj — 7T/128 and 6 = 7t/6. The signal is shown in 
Figure 3.16. For the time-dependent AR model while using the Fourier and Cheby- 
shev bases, best results are obtained with 33 and 26 basis functions respectively. 
The error variances in comparison with that found by non-parametric methods are 
given in Table 3.5. 

Table 3.5: Error variance in estimating the amplitude and frequency modulating 
signals. 

Modulating Non-Parametric Approach TAR Process 

Signal implemented with 


AM 

FM 


Non-linear Energy Hilbert Fourier Chebyshev 

Operator Transform Basis Basis 

1.69 x 1CT 3 1.41 x 10~ 3 3.41 x 10" 4 r 3.39 x 10" 4 

4.68 x 10 — 4 2.67 x 10“ 3 5.57 x 10~ 6 • 5.75 x 10~ 6 


The loci of the poles for the AFM process as estimated by the parametric method 
using the Chebyshev basis are shown in Figure 3.17. 

Since estimation of amplitude modulating signal is very sensitive to noise, the 
system is overdetermined and L = 7 is used. For estimation of the modulating 
signals a total of 32 Fourier bases and 24 Chebyshev bases are employed respectively. 
The performance of the proposed method in comparison with that of the non- 
parametric methods is shown in Table 3.6 and 3.7. 
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Figure 3.16: The AFM Signal y(n) = (1 4- kcosu a n) cos (9n 4- cos (w^i)), 
n — 1, . . . , 512; where, k = 0.5, u m = 0.50, u a = u>f = 7r/128, 6 = 7r/6. 
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Figure 3.17: The lod of the poles of the AFM signal. 



Table 3.6: Error variance in estimating the amplitude modulating signal of the 
AFM signal with overdetermined system. 
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Table 3.7: Error variance in estimating the frequency modulating signal of the 
AFM signal with overdetermined system. 
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3.4 Discussion 

Parametric modelling of non-stationary signals by the time-dependent ARMA 
(TARMA) model is investigated here. It is shown that the approach is generalised, 
and modelling by complex exponentials, AFM, etc. signals are special cases und«J 
the general framework. This fact is confirmed by both synthesis and analysis of 
different type of signals. In the general framework of TARMA model, the damped 
sinusoid, AM, FM and AFM models can be classified by the positions/loci of the 
poles of the process. In case of the damped sinusoidal model, the poles are located 
within the unit circle. When the damping factor is constant, the poles are located 
at fixed points, whereas for time-varying damping the loci of the poles are given by 
the radial lines. For the AM signal model, the loci of the poles are the radial lines 
intersecting the unit circle as shown in Figure 3.10 . For the FM signal model, the 
poles have only the circumferential movements on the unit circle as shown in Figure 
3.15. For the AFM signal model, both the radial and circumferential movements of 
the poles are observed and the loci of the poles intersect the unit circle as shown in 
Figure 3.17. 

Transient signals like nuclear magnetic resonance (NMR) [21] and electromag- 
netic pulse (EMP) [58, 45] responses are traditionally represented by the complex 
exponential models. Example 1 of section 3.3.2 (analysis) shows that an arbitrar- 
ily damped sinusoid can be represented more suitably by the time-dependent AR 
(TAR) model, with a damping factor as a function of time. The TAR model is 
found to be more flexible, and it gives a new way of accurate modelling of transient 
signals. 

As shown in the examples of section 3.3.2 (analysis), the parametric method of 
the TAR model clearly outperforms the non-parametric methods of estimating the 
modulating signals in the no-noise case. In case of a signal corrupted with additive 
white noise the parametric method excels in performance when suitable basis func- 
tions are used. The estimation of frequency modulating signals is found to be more 
robust than the estimation of amplitude modulating signals, when corrupted with 
additive white noise. This is because in case of amplitude modulation any error in 
estimating 0.2(71) at any instant ‘n’ gets accumulated as shown in (3.32) and (3.34). 

Conventionally, the speech like signals are modelled by the AM/FM signals for 
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the purpose of analysis/ synthesis [46 47 481 tv, 

13 [W, 47, 48 j. The presence of the AM and FM 

signals m speech production was reported earlier [281 Recentlv H, . 

i , , recently the experiments on 

speech resonance also indicate the presence of the AM and PM „• , „ 

I55 . 56] The PM signal is frequently encountered in sonar ap^LloT^ 

underwater environment, the rotating machinery (vie. submarine prop JJ et “ 

produces the AFM signal [54], These phenomena motivate us to develope th 

parametric methods for estimating the AM and PM signal with greater accuracy 

Another point to be mentioned here is that unlike the non-parametric methods the 

parametric approach can be extended to the case of multiple AM or PM sinusmds. 

Preliminary results indicate the potential of the parametric method and further 

work is needed to use these techniques in real life situations. 


It should be emphasised that in this thesis we have tried to find the r elation s 
among various existing non-stationary models and exploit the relations for better 
modelling and estimation. For the purpose of estimation of model parameters, how- 
ever, we have relied on the techniques which are already available. The shortcomings 
of the estimation procedures are discussed in section 5.2 in detail. 



Chapter 4 


ECG Signal Synthesis by 
Exponential AM Sinusoids 

4.1 Introduction 

The Electrocardiograph (ECG) is a non-stationary signal [2], The burst-like QRS 
feature attributes localised high frequency components in the ECG signal, mak- 
ing it distinctly nan-stationary in nature [60]. Though this feature of QRS wave 
has helped detection of the wave by filtering/linear prediction [29], it makes the 
modelling of the signal very difficult. 

Most of the work in modelling ECG are non-parametric in nature [61, 62]. It 
is well known that parametric modelling is superior to non-parametric modelling 
for improved accuracy and flexibility in analysis and synthesis. The improvement 
results because more information is provided about the signal [1]. 

Previous attempt to represent parametrically a segment of the ECG either by 
the impulse response of a pole-zero model or by Prony’s model was not very suc- 
cessful because the models had not considered explicitly the non-stationarity of the 
modelled signal or its pseudo-periodicity [63, 64]. This led to prohibitively large 
model orders, or more seriously, the model could represent only the local character- 
istics of the ECG signal. It has been, therefore, the main motivation of the present 
thesis to develop a parametric model for the ECG signal, which will represent its 
global nature. 

The time-dependent AR/ARMA model is the representative of the general class 
of non-stationary signals [2]. As such the model can be used for the ECG signal as 
well. However, the ECG has some distinctive features viz. i) its pseudo-periodicity; 
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u) different features of the constituent signals (P, QRS and T) representing actions 
of various parts of the heart [30]; etc. It will be insightful to know how the general 
time-dependent AR/ARMA model is restricted by the special features of the ECG- 
type signals. 

It is shown in this thesis that the amplitude modulated (AM) sinusoidal signal 
model is a special case of the time-dependent AR/ARMA model. The AM signal 
model can have the periodicity property, and can exhibit burst-like feature very 
well, when the modulating signal is an exponential function. Hence, it is proposed 
that one or more Ahd sinusoidal signal(s) can be employed to model separately 
each constituent signal of the ECG. The suitability of the developed model is then 
investigated for the ECG signal using analysis-by-synthesis te chni que 


4.2 Study of the Model 

4.2.1 Model Development 

In Section 3.2.2, it is shown that from a single order AR process (3.6) : 

*(„) = <' aX(n - 1] +W{U) i£n<n “ (4.1) 

k ad(n)x(n — 1) + w(n ) if n > no 

with a « 1.0 exp (J9), d(n) = exp (b sin a(n — no)) n > no and the noise variance 
approaching zero, one can get the sequence 


y(m) = Real [a; (no + m)] 


A exp (1 — cos am) cos (9m + (f>) 
\ a ) 


(4.2) 


where a; (no) = Aexp(j<f>) . The sequence y(m) is an amplitude modulated (AM) 
sinusoid with carrier frequency 9 and modulating factor given by exp (— ^ cos arnj . 
We term the signal y(m) as outcome of exponential amplitude modulation (EAM). 
The presence of exponential form helps in generating the burst-type periodic signals 
like the ECG. 


Analysis-by-Synthesis Rules 

The parameters of the signal model can be determined through the estimation of 
process parameters, as suggested in Reference [2]. A suitable set of basis functions 
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can be used to express the time-varying process parameters. We will, however, adopt 
here an analysis-by-synthesis approach for finding the parameters of the ECG signal 
directly. The approach will provide some insight regarding the composition of'the 
signal. 

In the analysis-by-synthesis approach, first note that the carrier signal fre- 
quency Q must be an integer multiple of the modulating frequency a to have similar 
waveshape in every period of the signal y(n). Though the ECG signal is only 
pseudo-periodic, we have considered the ideal case where the modelled signal be- 
comes a periodic extension of one period of the ECG signal. It is done to make the 
synthesis procedure simple. Otherwise, making 9 and a time-varying while keeping 
their harmonic relationship intact, the actual ECG data can also be modelled. 

Under the condition of strict periodicity for y(m), (4.1) can be rewritten as , 


y(m) = A exp 


(1 — cos am) ^ cos ( kam + <j > ) 


(4.3) 


where k is an integer. 

To generate a complicated signal pattern more than one AM sinusoids are 
needed. The most prominent feature of the signal is first fitted by one AM wave, and 
once fitted, this part of the original signal is removed. We continue this process on 
the residual signal till all the features of the signal under investigation is exhausted. 
For the purpose of fitting an EAM sinusoid, the following steps are employed : 


1. We start with some initial guess of A , b and k. At the starting (j> — 0, and 
a = ~ , where T is the period of the signal under investigation. 

2. The value of the parameters 6 and k are to be altered to get proper shape 
of the wave. The phase is also changed to decide the sign of the peak and 
skewness of the hump. So, the parameters 6, k and <j> are changed iteratively 
till the shape of the synthesised signal comes close to the target feature of the 
ECG signal e.g. the QRS complex. 

3. Then, the signal is scaled by a factor c and the origin is shifted to mo, so that 
the target feature of the original signal coincides with the synthesised one. If 
the fitting is not satisfactory at this stage, we have to go back to step 2 and 
retry. The synthesised signal z(m) is expressed as, 


z(m) = c y(m + mo) 


(4.4) 
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4. Once the synthesised signal is ready, it is subtracted from the signal under 
investigation. In this way, we get the residual signal. If the residual signal 
still has some prominent feature, we have to go to step 1 to start with a new 

EAM sinusoid. When all features of the signal are exhausted, the process is 
stopped. 

4.2.2 ECG Signal Fitting 



Figure 4.1: The original ECG signal, Lead II configuration, 180 samples sampled 
at 250 Hz (solid), and synthesised ECG signal (dash). 

The Lead II ECG signal of a normal male patient collected for one period and 
sampled at 250 Hz is shown in Figure 4.1. The T and QRS waves are synthesised 
by two separate AM sinusoids. For P wave, two more AM sinusoids are needed, and 
the corresponding waves are named as PI and P2. A d.c. signal is also needed to 
match the d.c. bias of the original ECG signal. The successive waves are matched 
in the following order : T, QRS, PI, and then P2. These waves are shown in Figures 
4.2 - 4.3. While fitting the T wave, the d.c. bias of the signal is adjusted. 

To synthesise the ECG signal components, we start with |a| = 0.9999, cr w = 10 28 
and no = 15000, as defined in (4.1) . The d.c. bias is found to be 500. The values 
of the other parameters for the constituent waves T, QRS, PI and P2 are given in 
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the following table. 


Table 4.1: Parameter values for the ECG signal synthesis 


Parameters 

T wave 

QRS wave 

Pi wave 

P2 wave 

a 

0.03808 c 

0.03808 c 

0.03808 c 

0.03808 c 

b 

0.19 

1.1 

0.1 

0.1 

k 

2 

7 

1 

1 

4> 

0.25tt c 

1.157t c 

1.67t c 

7T C 

, mo 

5 

54 

155 

85 

c 

0.78557 

4.778 x 10 -21 

40.0 

15.0 

A 

7.258 x 10~ 2 

6.917 x 10~ 2 

2.505 x 10“ 2 

2.505 x 10 


The synthesised ECG signal is shown in Fig 4.1, which closely resembles the 
original signal, over the period of consideration. As the synthesis technique has 
produced a periodic extension of one period of the ECG signal, we compute SNR 
and correlation coefficient over the period concerned. The SNR is calculated to be 
9.06dB, when SNR is defined as, 


SNR = 10.0 log 


E"., (*„(i)) 2 

S.1 (WO - WO ) 2 


where x ors (i) and Xsy n (i) are the original and synthesised ECG signal at i-th instant, 
and n is the period of the ECG signal ( = 165 samples ). The correlation coeffi- 
cient (p) between the synthesised and original ECG is calculated to be 0.933. The 
correlation coefficient is defined by 

— m org)(. :I 'syn{' n ') m syn) 

^/S”=l( x ws( n ) — m org) 2 T>t=l( X *yn{ n ) ~ m syn)' 1 

where n, x^i) and x* yn (i) are as defined earlier, and m^ and are the mean 
of original and synthetic ECG signal respectively. The generating process can be 
represented as shown in Fig 4.4. Equivalently, Fig 4.5 provides the complete process 
representation of the signal. Note that the functional block, Real[ • ] does not appear 
in the latter diagram. 
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Figure 4.2: The synthesised QRS wave (solid), and synthesised T wave (dash). 
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Figure 4.3: The synthesised components of P wave, PI (solid) and P2 (dash) 
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Initial condition* 
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Figure 4.4: The block diagram of the ECG signal synthesiser by first order complex 
AE models. 
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Figure 4.5: The block diagram of the ECG signal synthesiser by second order real 
AR, models. 
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4.3 Remarks 

The method presented here provides a novel way to parametrically model the ECG 
signal, and also its constituent waves. It is known that the P and QRS waves 
are caused by the electric potentials generated as the depolarisation prior to the 
atria and ventricles respectively. Whereas, the T wave in the ECG is caused by 
the potential generated as the repolarisation of the ventricles. Therefore, it is 
meaningful to model the constituent waves separately as it is done in the present 
approach, because they correspond to the separate actions of the various parts of 
the heart [30]. 

It is apparent that the ECG generating model as shown in Figure 4.4 or 4.5 
gives an equivalent process representation of various functions of the heart. As 
almost all the serious abnormalities of the heart can be detected by analysing the 
ECG wave from different leads, this model can be expected to provide insight in the 
functioning of the heart, and may be helpful for diagnosis [30]. It is clear however 
that considerable amount of further research will be needed in this direction before 
the model can be employed for clinical diagnosis. 

Until now, the data compressions for the ECG are implemented mostly by non- 
parametric approaches because the parametric methods, in most cases, did not 
preserve the waveform [65]. The method presented here being a synthesis technique, 
guarantees the preservation of waveform. Moreover, the use of a small number of 
parameters, (only 25 parameters required) promises the application of the model 
for superior data compression, but it can only be asserted after further research. 

For the analysis of the signal, direct evaluation of model parameters leads to a 
set of highly nonlinear equations. As an alternative, one may estimate the equiva- 
lent process parameters by the method as shown in Appendix A. The equivalence 
between the two approaches as shown in this thesis is an important concept. The 
concept can be conveniently utilised for parameter estimation of the signal model. 

The analysis-by-synthesis approach as demonstrated here is simple to imple- 
ment. The approach clearly provides good insight about modelling of the ECG 
signal [30]. 


Chapter 5 


Conclusion &; Scope for Future 
Work 


5.1 Conclusion 

In the study of TARTG model, we have found that both the time-varying AR co- 
efficients and the time-varying gain are important in their own ways. The first 
set represents the spectral waveshape variation over time and the latter keeps the 
information of time variation of the spectral gain. Obviously, none can represent 
the information of the. other. The estimation of gain allows us to obtain a faithful 
estimation of the evolutionary spectrum, and in this way, the interpretation of the 
evolutionary spectrum as the energy density in the time-frequency plane becomes 
apparent. 

In parametric modelling, the parsimony of the model is an important issue. If 
a model requires a large number of parameters, then it is clear that the model is 
not properly fitting into the signal [66]. In this case, search should be continued 
for more suitable model. In modelling the speech phonemes, we have found that 
the QSAR and TAR models require a large number of parameters compared to the 
ARTG model, while the signal-to-prediction error ratio varies only slightly. So, the 
ARTG is a more suitable model for speech phonemes, which can be directly used 
in speech synthesis. In case of speech coding, often the QSAR model is used, and 
thereafter, the redundancy of parameters are removed by different vector quanti- 
sation techniques [35]. If in continuous speech, the phoneme boundaries or more 
precisely, segmentation of frames with different spectral waveshapes can be done, 
the ARTG model can be used for speech coding. The alternative model can save 
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a lot of computation for searching of the code-book. In the same way, the ARTG 
with proper segmentation will be a better choice for linear prediction (LP) based 
speech recognition algorithm. 

Even in case of the TAR model (with constant gain) when used for coding, the 
estimation of gain should be used for normalisation. This can be illustrated by 
an example as follows. Let us consider two different stationary AR processes of 
different signal strengths placed one after another. When a TAR model is fitted for 
the combined sequence using a finite number of basis functions, the obvious result 
of estimation of parameters by minimising the unweighted error function will be a 
close fit in the portion where the signal strength is high, and will not follow the 
portion of the signal where the signal strength is low. As in a TAR model any frame 
under consideration can have vowels and consonants, the effect will be s imil ar if the 
signal is not normalised using the estimated time- varying gain. 

The investigation on the movement of poles of the time-dependent RTF model 
has revealed that the damped sinusoidal model, AM sinusoidal model, FM sinusoidal 
model or AM/FM sinusoidal model are only special cases of the generalised model. 
Each of these models can be generated by putting proper restrictions on the poles 
of the time-dependent RTF model. This knowledge gives us a unified procedure for 
estimation of parameters. This leads to more flexibility/accuracy of modelling of 
damped sinusoids and better accuracy in estimating the modulating signals in case 
of modulation models. 

The last part of the thesis is devoted to modelling of the ECG signal. The 
synthesis technique presented has given insight about the ECG modelling. The 
developed model is promising and it may be useful in explaining various func- 
tions/malfunctions of the heart. 

5.2 Scope for Future Work 

The study on parametric modelling of non-st ationary signals as shown in this thesis 
has shown great potential for various applications in future. Before harnessing the 
advantage, however, a few hurdles have to be crossed. In the present work of unifi- 
cation of various linear non-stationary models and classification of non-stationary 
signals, the ‘standard’ method of estimation through basis functions is used for each 
AR coefficient of the TAR model. The method, however, has several drawbacks and 
some of those are demonstrated in section 3.3.2 (Alalysis). Therefore, one has to 
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find suitable, solutions to a number of issues before the modelling technique can be 
applied in real-life solutions. The issues are : 

• Choice of the basis function 

• Number of the basis functions to be used 

• Order of the model 

Although these issues are addrssed separately, these are highly dependent on one 
another. The choice of basis function refers to what should be the natural choice 
of basis. For example, when the time-variation of the parameter is linear, the 
sinusoidal basis is quite inappropriate and its choice may lead to a very large basis 
order. So, before employing any set of basis function, one must have an idea about 
the time-variation pattern of the parameter under consideration. The issue of 
optimal basis order selection arises next when a set of basis functions is selected. 
The order of the model and the number of basis functions should be employed are 
also inter-dependent. As the error in fitting the model (i.e. prediction error in 
case of the AR model) can be minimised either by increasing the number of bases 
or by increasing the model-order, the optimal choice of these quantities are very 
difficult. The question may be complicated for a multi— component signal where 
only one component is prominent at a time, e.g. the ECG. Apperently, we can have 
a large model-order with a small number of basis functions or a small model-order 
with a large number of basis functions. So, the conventional root-mean-square error 
judgement criterion can not be expected to give a prudent answer to this question. 
For this purpose, some new approach and new theory should be developed for non- 
stationary signals. This should be done while utilising the information about the 
production/perception mechanism of the signal in concern. Once these questions 
are reasonably answered, attention can be diverted to make the algorithms more 
robust from noise. This remains to be the motivation for future research and the 
suggestion for the future researchers in this field. 



Appendix A 


Estimation of Time-Dependent 
AR Coefficients 


The time-dependent autoregressive (TAR) model is represented by the equation, 

x(n) + ai(n)x(n — 1) 4 b Op(n)x(n — p) = w(n ) (A.l) 

for n = 1,2 where p is the autoregressive model order, w(n) is the white 
noise sequence and, a^n) is the i-th AR coefficient at an instant ‘n’. 

The coefficient a.j(n) is assumed to be ‘smooth’ in the sense that if the first 
derivative of the coefficient may be arbitrarily great, the higher order derivatives 
necessarily vanish. So, the coefficients a;(n) can be approximated by a set of basis 
functions. The time-dependent coefficient aj(n) is expressed as, 

m 

a >( n ) ( A - 2 ) 

;=0 

where fj (n) is the value of j th basis function at an instant ‘n\ Let us consider the 
quantity a.i(n)x(n — i ) which can be written as, 

( ™ \ 

ai(n)x(n — i) = 22aijfj{n) x(n-i ) (A.3) 

Vi=o J 

m 

= 5Za ii (/ i (n)as(n-i)) 

j= o 

mT 

Let us define the vector X(n,i ) = f 0 (n)x(n — i), ..., f m (n)x(n — i) ~ , where 
symbol e T” stands for transpose of a vector. With this definition (A.3) becomes 
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Using (A 4), the TAR process (A.l) pplp be written as 

x ( n ) + ^ X T (n, 1) • • • X T (n,p) 9 = w(n ) (A.5) 

where 

a — ' ~ T 

^ O'lO ‘ ‘ * 0,\ m 0.20 " * ' &2m * ’ ’ npO * * * &pm (A.6) 

The interpretation of (A.5) is immediate. The parameter vector 9 is the vector 
of regression of the non-stationary process x(n) on the past samples of the vector 
X(n,i ) for i = lj . . . ,p. Here lies the advantage of basis parameterisation : a 
linear time varient or non-stationary problem becomes a linear time-invariant or 
stationary problem by replacing a scalar process by a vector one. 

The innovation process w(n ) of the TAR can also be viewed as a prediction error 

w(n) = x(n) — x(n) (A.7) 

where x(n) = E(x(n)\x(n — 1), • • • , x(n—p )) is the prediction of x(n) given the past 
values x (n — 1), ... , x(n — p), i.e. 

x(n) = - ^ X T (n, 1) • • • X T (n,p) 9 . (A.8) 

It is therefore me ani ngful to minimise the variance of this error w(n) and hence 
we get [2], 

E [(grad ff u;(n)).w(n)] = 0 (A.9) 

The gradient can be easily computed as 

grad e to(n) = X T (n, 1) X T (n,p) ^ (A.10) 

We therefore obtain the op timum vector 9 as the solution of the Yule-Walker type 
equation, 

(~X(n, 1)" > > 

E i " X T (n, 1) ... X T (n,p)] 9 = -E : x(n) 

{ X{n,p) ) V. X(n,p) _ / 

V " ' (A. 11) 

Remark 1: If the coefficients are constants i.e. m — 0 and fo(n) = 1, then (A.11) 
becomes the Yule-Walker equation of AR(p) process. 
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Remark 2: If the process x(n) is purely auto-regressive but, w(n ) has a gain g(n), 
it is no longer significant to minimise the prediction error.- However, (All) 
remains valid in the sence it expresses nothing but the orthogonality of the 
innovation w(n) and this orthogonality still makes sense, even if g(n) is not 
constant [2]. 

As (All) is obtained through minimisation of the variance w(n), which is char- 
acterised as a stationary and white noise process, it is possible to use in (All) 
the ergodic estimator for expectation, even if the original process x(n) is non- 
stationary. This allows us to replace the expectation by a summation over a time 
interval r = [p + 1, N] [2]. 

When the signal is embedded in additive white noise, an overdetermination in 
the matrix equation (All) can give better estimate [32]. In this case, (All) takes 
the form 

X(n,l) ' \ ^ "" X(n, 1) 1 \ 

E X{n,p) ' X T {n, 1) ••• X T (n,p) ' 6 = -E X(n,p ) x(n) I 

^ L X(n,L) ) ^X(n,L)_ j 

(A. 12) 

where L > p. The extra set of equations results from the fact that the innovation 
w(n) at an instant ‘n’ is also orthogonal to the past samples X(n,p+1), . . . ,X(n,L). 
In this case, the time interval r over which the summation replaces the expectation 
becomes [L -f 1, N}. 

For the purpose of representing the time-variation of the AR coefficients, the 
Fourier and Chebyshev basis functions are employed in this thesis. The k-th Fourier 
basis is defined as 

cos for even value of k, 

/*(«•) ^ sin for odd value of k, 

where 1 < n < N. The k-th Chebyshev basis is defined as 

f k (n ) = cos (fcarccos (Z — 1)), 

where Z = and 1 < n < N. 

— I — — 



Appendix B 

Fourier Cosine Series with 
Half— Range Expansion 


Let a continuous time signal f(t) be defined on an interval 0 <t <T. Then Fourier 
cosine series with half range expression is of the form [67, Chapter 10] 


00 / \ 

fit) = Co + Y, COS (jY *) » 0 < t < T 

(B.l) 

and the coefficients are given by 


co = 

Cn = ~ jT /(*)cos(^t)dt,n = l,2,.... 

(B.2) 

(B-3) 

For the discrete time signal /(fc), k = 1, the sampled version of the con- 

tinuous time signal f(t ) for 0 < t < T, the Fourier cosine series 
expansion of f(k ) is defined as 

with half range 

/(A:) = co + £]c„cos k = l,2,...,N. 

n=l v 

(B.4) 

and the coefficients are expressed as 


Co = TrYlfW 

Jt=l 

(B.5) 

* = ^|;/Wco S (^),n= 1 ,2 : W. 

^ k—1 

(B.6) 
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Appendix C 


Representation of a Sinusoid by 
the Stationary AR Model 


Let us consider a real sinusoid x(n) of the form, 

x(n) = A cos (27r/ 0 n + (f> ) (C.l) 

Strictly speaking the use of any time series model for a sinusoid is inappropriate. 
The AR(2) model is given by 

x(n) = — aix(n — 1) — a, 2 x(n — 2) + w(n) (C.2) 

where w(n ) is the input white noise with variance a~. This model can not generate 
a sinusoidal time series. A recursive difference equation that can generate (C.l) for 
N > 0 is 

x(n) = 2 cos (27r/ 0 )x(n— 1)— x(n-2)+Acos^6(n)-Acos(27r/ 0 — (j>)6(n— 1) (C.3) 

where x(— 2) = x(— 1) = 0 and S(n) is the ‘dirac delta ’ function. The poles of the 
all pole filter (C.3) are the solution of the equation 

A(z) = 1 — 2 cos (27rfo)z~ 1 + z~ 2 = 0 


or, 

Zp = exp (±j2irf 0 ) (C.4) 

and hence (C.3) represents an unstable filter. 

Inspite of these difficulties, it is still possible to model the auto-correlation func- 
tion (ACF) of (C.l) as the limiting form ACF of a real AR(2) process. The ACF 
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of x(n) as given in (C.l) is 

A 2 

T xx(k) = — cos (2-Tr/ofc) (C.5) 

This ACF can be approximated by ACF of an AR,(2) process with poles at 
r exp (±27 t/ 0 ) if r —» 1 <r 2 -» 0 in such a way that r„( 0) = 4r- 
The ACF of an AR(2) process can be shown as 



In this way a sinusoid may be modelled by an AR(2) process, at least as fax the 
ACF is concerned. For M real sinusoids a real AR(2M) model is required, and for 
M complex sinusoid complex AR(M) is needed [1, Chapter 5]. 



Appendix D 
A Long Proof 


m-hriQ 

Y sin a(i — n 0 ) 

i=l +no 


m 

Y sin ai 

»=i 


= Y, s i n a i 

»=o 

m 

= Im Y ex P (j&i) 

4_i=0 

_ "exp {ja(m + 1)) - 1" 
exp (joe) - 1 

cos a(m + 1) — 1 + j sin a(m + 1) 

= im . . 

cos a — 1 + j sin a 

__ ^ "(cosa(m + 1) — 1 + j sin a(m + 1)) (cos a — 1 — jsina)' 

(cos a — l) 2 + sin 2 a 

cos a sin a(m + 1) — sin a(m + 1) — cos a(m + 1) sin a + sin a 
cos 2 a + 1 — 2 cos a — sin 2 a 

[sina(m + 1) cos a — cosa(m + 1) sin a] + [sin a — sin a(m + 1' 

2 — 2 cos a 

sin am + 2 cos (a + 2 ?p) cos ^ 

4 sin 2 § 

cos sin zf- + cos (a + s y)^ 

= 2 sin 2 1 

When a is small, the above expression can be reduced to a more simple expres- 
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sion by substituting summation by integration. The derivation is given below 

m 

Esin (ad) = ^sin(ai) 
i=1 »=o 

rm 

- J o sin (at) dt 


- ~[- cos (^)3" 

1 

[1 — cos (am)] 


a 
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